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Abstract 

The Liouville-space Green function formalism is used to compute the current density profile 
across a single molecule attached to electrodes. Time ordering is maintained in real, physical, 
time, avoiding the use of artificial time loops and backward propagations. Closed expressions for 
molecular currents, which only require DFT calculations for the isolated molecule, are derived to 
fourth order in the molecule/electrode coupling. 
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I. INTRODUCTION 



There is a considerable interest in measuring 1-14 and computing^ " 18 i 19 ~— currents in 
single molecules attached to electrodes. The current- voltage (J — V) characteristics of such 
devices are strongly nonlinear, and should be recast in terms of high order conductivi- 
.j-jggi^isj^ Controlling molecular currents has important applications to nanodevices such 
as molecular wires and rectifier o 12 ! 13 ' 14 . 

The first-principles computation of molecular currents poses a serious challenge that was 
addressed by many approache d 17 ' 18 . Standard quantum chemistry packages are not designed 
to compute current-carrying states. Non-equilibrium Green function techniques widely used 
for computing currents in macroscopic systems, e.g. semiconductors^ 2 ^, were adopted for 
single-molecule current o 1 ^ 1 ^ 17 . The underlying diagrammatic perturbative expansion re- 
quires the computation of extra information in the intermediate steps, which is not needed 
for computing the current. This includes the dependence of the Green function on two times, 
involving coherence between states with a different number of electrons^ 7 -. In a different ap- 
proach the time dependent charge distribution in a donor-acceptor junction was studied 
by calculating the tunnelling currents at the Hartree-Fock level in terms of the donor and 
acceptor states^. 

An ab initio algorithm for calculating molecular currents developed by Lang 21,22 is based 
on the self-consistent solution of a system of four equations: (1) The Kohn-Sham integral 
equation for the wavefunction of the molecule + electrodes (Eq. (3.1) in Ref.—), (2) Equation 
for the potential difference between the complete (molecule + electrodes) system and the 
bare electrodes (Eq. (3.2) in Ref.—), (3) Equation for the Green function of the molecule 
in an effective potential given by the sum of the exchange-correlation and the electrostatic 
potential (Eq. (3.13) in Ref.—). The latter is determined by (4) the Poisson equation with 
boundary conditions setting the potential to be equal to the chemical potentials difference 
of the electrodes (Eq. (jHJ) in this paper). These equations, together with the boundary 
conditions, yield the wavefunction of the entire (molecule + electrodes) system which can 
be used to compute the finite current. This algorithm was applied for calculating the current 
in a chain of three Al atoms with the substitution of one atom by sulfur, and for a chain of 
ten carbon atoms. The electrodes were modelled as a semi-infinite ideal metal (homogeneous 
electron gas or Jellium model). The same approach was subsequently applied to calculate the 
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current in a benzene-1, 4-dithiol molecule connected to ideal metal electrodes, and the effect 
of inserting a single gold atom between the molecule and a metallic surface was investigated^. 

Mujica, Kemp, Roitberg and Ratner— proposed to solve self-consistently the Hartree- 
Fock equation for the Green function of the complete (molecule + electrodes) system and 
the Poisson equation for the electrostatic potential, with the same boundary conditions of 
Lan g 21 i 22 . The molecule was described by the Hubbard Hamiltonian, the electrodes were 
taken to be an ideal metal, and the current is computed from the Green function. This 
approach is easier to implement compared to the density functional theory (DFT) 21 i 22 , since 
it only requires the solution of two rather than four self-consistent equations. However, 
DF T 21 i 22 takes into account exchange-correlation effects more rigorously. 

Kosov applied a variational Kohn-Sham density-functional (DFT) equation for electrons 
in molecular wires (benzene-1, 4-dithiolate molecule covalently bonded to two Gold elec- 
trodes)^. By variational minimization of the total energy, where an arbitrary given current 
is imposed as a constraint, a ten atom Gold cluster participates in the bonding with a 
molecule. The device contains 32 atoms. This results in a system of two self-consistent 
equations, similar to the Hartree-Fock approach 2 ^, which is simpler than the equations of 
Ref. 21,22 , yet treats exchange-correlation effects as rig orously as DFT^ 2 -. 

All of these approaches, which employ DFT for the molecule + electrodes, require the 
solution of a system of self-consistent equations, but do not give closed expressions for the 
current. 

In this paper we compute the current profile in a molecule perturbatively in the ratio 
of the coupling of the molecular electrons with the electrodes and the separation between 
Kohn-Sham orbital energies. Our approach yields a closed expression for the current density 
(|j(r)^ in terms of Kohn-Sham orbitals of the molecule alone, which enter various Green 
functions in the frequency domain. The molecule/electrodes coupling potential of Lang's ap- 
proach 2 ^ is treated to the lowest order required to yield a finite current. We use a Liouville 
space (superoperator) many-body Green function technique, which provides a convenient 
description of non-equilibrium effect a 2 ^ 2 ? . The Liouville space TDDFT (Time-Dependent 
Density Functional Theory) formulation of non-linear response based on the single electron 
density matrix was developed in Ref . 28 i 29 , and subsequently extended to superconductors^. 
One advantage of working in the higher dimensional Liouville space is that we need only con- 
sider time ordered quantities in real (physical) time and Wick's theorem therefore assumes 
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a particularly compact formal; in contrast to Hilbert space Green function perturbation 
theor y 24 i 31 i 32 no special contours or analytic continuations are necessary. 

II. THE MODEL HAMILTON! AN 

We shall partition the electronic Hamiltonian of a molecule connected to two electrodes 
A and B, in the form 

H = H + V, (1) 

where Ho represents the noninteracting molecule and electrodes, and V is the 
molecule/electrodes coupling 

V = [dx [ dx'U{x,x') W{x){^ A {x') +ip B {x')) + {ip A] {x') +i) B] {x'))i){x)] . (2) 



The spatial coordinates x and x' run over the molecule and the electrode regions, respectively 
(x = (r, t) is a coordinate of electrons within the molecule, and x' = (r', t) is a coordinate of 
electrons in the electrodes region). U(x,x') is the coupling potential between the molecule 
and electrodes; ip{x) is the field operator of electrons in the molecule, whereas -ip A (x') and 
i/j b (x') are the field operators of electrons in the two electrodes. These operators satisfy 
Fermi anti-commutation relation o 33 ! 34 

^(ri)^ t (r 2 ) + ^(r 2 )ip(ri) = 5(r 1 - r 2 ); 
^(ri)^(r 2 ) +^(r 2 )^(ri) = 0; 

^ t (r 1 )^( r2 ) + ^t( r2 )^t( ri ) = . (3) 

We assume that the electrodes interact with the molecule only at two points of contact 
(r^ and Tq ). Kosov considered the coupling between the molecule and the electrodes rep- 
resented by 10 Gold atoms^; participation of only a restricted domain of the electrodes in 
the molecule/electrode coupling was also assumed in the "extended molecule" model^. We 
then have 

f/(r,r') = 0(r)[ ( 5(r'-r^) + 5(r'-^)], (4) 

where 0(r) is defined as the potential, induced by the difference between the charge density 
of the isolated molecule £o( r ) an d of the molecule coupled to the electrodes ^(r) 1 -. 0(r) is 
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obtained by solving the Poisson equation in the molecule (r A < r < r B ) (Eq. (5.3) in Ref.— ) 

V 2 0(r) = -e(g(r) - g (r)), (5) 

with the boundary conditions 0(r^) = fiA and 0(i"q ) = /ie. Here ha and /^s are the chemical 
potentials of electrodes A and B, respectively, W = fiA — is the external voltage, and e 
is electron charge. 

The Poisson equation for the total electrostatic potential in the molecule coupled to the 
electrodes U e i(r) i a 17 i 19 

V 2 U el (v) = -eg(r), (6) 

with the boundary conditions f/ e ;(r^-) = \ia and U e i(r^ ) = \i B . In act the boundary condi- 
tions should reflect the change of the electrostatic potential in the interior of the left (right) 
electrodes when current is flowing, rather than the chemical potential difference of electrodes 
(Eq. (5.4) in Ref.-^). The latter between the last two is, in general, smaller than the applied 
voltage^. This difference should be negligible for large electrode spacings, for which the 
current is small, but should be noticeable for small spacings. Our boundary conditions thus 
hold in the small current limit (large electrode spacings). 

To lowest order in the molecule/electrodes coupling we use the unperturbed charge density 
of the isolated molecule in the r.h.s. of the Poisson equation ^(r) = Qo(r) , and Eq. © gives 

V 2 ^(r) = -e5> Q (rK(r), (7) 

a 

where the sum runs over all occupied molecular Kohn-Sham orbitals </2 Q (r). Using g(r) = 
q (t), we obtain from Eq. © (Eq. (5.6) in Refill) V 2 0(r) = 0, 0(r£) = fi A , 0(r£) = \i B . 
By solving this equation, we obtain that the coupling potential is linear in r 

0(r) = (r£ - r^)- 1 (/^ - ^ + {^a ~ Hb)t) . (8) 

There are two contributions to the charge density in the Poisson equation (Eq. (jB|)): The 
charge density of the isolated molecule and the variation of the charge density with the 
electrodes due to the coupling induced by the external voltag e 1 ^ 21 . To lowest order in the 
coupling potential we neglect the latter. And the electrostatic potential in the molecule is 
calculated using the unperturbed charge density (Eq. (JJJ)). We thus neglect all screening 
effects and the potential profile across the junction is assumed to be linear (Eq. (jSJ)). This 
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holds for low voltages. The electronic structure of the molecule enters in H . The external 
voltage, which is the difference between the chemical potentials of the electrodes (W = 
Ha — ^b), enters in the Hamiltonian of the molecule/electrodes coupling V in Eq. (JTJ through 
the coupling potential U (Eqs. (jH) and (jHJ)). Our perturbative approach holds for low 
voltages, where provided the coupling potential of the molecular electrons with the electrodes 
(the external voltage (Eq. (0))) is much smaller than the difference between Kohn-Sham 
orbital energies. 

It corresponds to lowest order expansion in the molecule/electrodes coupling potential of 
Lang's approach 21 . 



III. PERTURBATIVE CALCULATION OF THE CURRENT 



Our goal is to compute the expectation value of the current operator 



ieh 



J(x) = - ^(V r -V,)^(x)^(x] 
zm 



(9) 



where e and m are the electron charge and mass, respectively. 

We start with the interaction picture expression for the expectation value of J(x), using 
the partitioning (Eq. of the Hamiltonian 27 . We define the Liouville operators L = 
L + V- corresponding to Eq. ((TJ) where L = (-Ho)- i-e., L X = H X — XH . The zero 
order time evolution operator is defined as 



Go(i~2,n) = 0(t2 - Ti)exp 



--L (t 2 - n) 
n 



(10) 



where 9(t) is the Heavyside step function. Throughout this paper we denote operators in 
the interaction picture by a tilde (~), i.e. 



= Gl(r, 0)A u g (r, 0) 



(11) 



-, or v = L, R. 



The current is given by 



\T J + {x)exp 



n 



drV. 



T 



(12) 



where the Liouville space superoperators with indices + and — are defined in Appendix IXF 7 . 
and the time-ordering operator in Liouville space T is defined in (Eq. (jBll|) ) 27 . This natural 



time ordering chronologically follows the various interactions with the electrodes. The precise 
order in which superoperators appear next to a T operator is immaterial since at the end 
the order will be fixed by T. T before an exponent means that each term in the Taylor 
expansion of this exponent should be time-ordered. 

The current can be also expressed in terms of the molecular electron density matrix 
p(r, r', t) in the interaction picture 



p(r,r, t)=Tr ip(r, t)${r, t)p 



(13) 



where po is unperturbed many-electron density matrix of the isolated molecule uncoupled 
to the electrodes. The expression for the current is 



J(x)) = -^(V r -Vr)p(r,f,t) 



(14) 



A finite current requires at least two interactions with each electrode. We, therefore, 
computed the current to fourth order in V 

ieh 



where 



S(x, x) 



J{x) 



i a 



2m 



(y f -v T )s(x,x) 



(15) 



4! U 



dr\ I dr 2 / dr 3 / rfr 4 

3 J — OO J — OO J — OO 

dr\ J dr'i J dr 2 J dr'2 J dr% J dr'3 J dr 4 J dv\ 

fp(x, x, t)VL(r 4 , r4)y_(r 3 , r 3 )y_(r 2 , r 2 )VL(r 1 , n) 



(16) 



where p(x,x,t) = ip(x)ip Ji (x). 



We next expand S(x, x) in Eq. (fT6|) in "left" and "right" operators in Liouville spac 



:e2£. 



To that end we introduce the field operator corresponding to both electrodes 



(17) 



S(x, x) 



drA 



Substituting Eq. into Eq. (fTTj) . we get 

J dr\ J dr\ J dr 2 J dr'2 J dr% J dr'3 J dr 4 J dr\ 

U (x 4 , x' 4 )U(x 3 , x' 3 )U(x 2 , x' 2 )U{x u x' x )I{x u x[, x 2 , x' 2 , x 3 , x 3 , x 4 , x 4 , x, x), (18) 



dr\ I dr 2 / dr^ / ui 1 

00 J —oo J — oo J — OO 
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where I is given by Eq. (|A4|) . In Appendix El we use superoperator algebra to bring it 
to the form of Eq. (jA8|) . Factorizing that equation, using Wick's theorem for fermion 
superoperators (Eq. (74) in Ref.— ), we obtain 



^- v v ,u K =0,l; i,j =0,1 



rj,K=0,...,4; k,l=0,l 

^,0^ K ,l4,l)( 1 _ ^«,0^ K ,l4,l)] , (19) 



where 5 is the Kronecker symbol: for u v = ^ = ?/>l, and for u v = 1 ^ = t/jr; for k = 1 
ip Mk (x k ) = ip M (x'), for A; = ip Mk (x k ) = if)(x); for z = ip ai = if), and for i — 1 *0 ai = ift . 
P = 1/2 fn^ + n*^) is a number of sum of permutations of ipL and tpu, where n d ^ is a 
number of contractions in the r.h.s. of Eq. (|19|) (T (if)L(x n )if>L(x m ))) , when n ^ m, and n# 
is a number of contractions in the r.h.s. of Eq. (fTT?|) (T{ip R (x n )ipn(x m ))) , when n ^ m. 

We assume that the molecule contacts each electrode at a single point (r^ and Tq ). The 
coupling potential between the molecule and electrode is 

f/(r,r')=E^<(r)^(r'), (20) 

a/3 

where U(r, r') is determined by Eq. y? Q (r) is the a'th Kohn-Sham orbital of the molecule, 
and 

Uap = JdvJ dv'U{v, r')^(r)^(r') = (^(r£) + <p* (r»)) / rfr0(r)^ Q (r), (21) 

where 0(r) is determined by Eq. (jHJ). Below we denote U a = £/ aQ . 

Substituting the expansion of the frequency domain Green functions and coupling po- 
tential in Kohn-Sham orbitals (Eqs. (jC4|) and (J2U)0 in the expressions for I ( Appendices IE1 
and lC]) . and substituting / in Eq. (fTHj) and 5* in Eq. (fTH|) . we obtain 

w n\ e i f duo f duJi f du 2 r , \ 

(22) 



mS 3 v " r " r i 2vr J 2vr 7 2tt 

JV 21 



£S^( r )^( r ) 

/3=1 i=l 

Eq. ()22)1 is our final expression for the current. The twenty one terms ij contributing to the 
current are given in Appendix ID1 

If the external voltage is zero, \ia = = Ho- Setting /i = 0, Eq. (jHJ) gives 0(r) = and 
U (r, r') = everywhere in the molecule. All twenty one terms Jj contributing to the current 
Eqs. flm)- (ID221) then vanish. 
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IV. DISCUSSION 



In this paper we have derived a closed expression (Eq. (|2*2j) together with Appendix|Dj) for 
the current profile within a molecule attached to electrodes, in terms of Green functions of 
the isolated molecule. The external voltage W = ha — Hb is equal to the difference between 
the chemical potentials of electrodes A and B. Only one characteristic of the electrodes, 
the Fermi energy levels fx a and yU#, enters their Green functions, the molecule/electrodes 
coupling potential and the final expression for the molecular current (Eqs. (jG15|) - (|C19|) ). 
This expression only requires DFT calculations for the molecule alone (which enters in the 
Green functions of the molecule in the expression for the current) rather than the molecule 
and electrodes as used in Refs.— ~— . The derivation employs the Liouville space Green 
function technique^, and the current is expanded to fourth order in the molecule/electrodes 
coupling our formula for the current (Eq. (|22j)) should agree with Lang's approach 21 to 
lowest order in the molecule/electrodes coupling potential. The closed expressions for the 
current obtained in the paper, compared to the previous considerations, is not the result of 
the new methodology, but stems from the perturbative considerations. 

We have established a one-to-one correspondence between Liouville space and Hilbert 
space Green functions for fermions (Appendix El and AppendixlHJ). For completeness, similar 
relations for Boson Green functions are given in Appendix |E| The only difference is the 
definition of the time-ordering in Liouville space (Eq. (jBlljl for fermions vs. Eq. (|E1|) for 
bosons). 

Liouville space provides a fully time ordered description of the current since we only need 
to propagate the density matrix forward in time. In contrast, the wavefunction in Hilbert 
space involves both forward and backward propagations. The choice is between following 
the ket only, moving forward and backward or following the joint dynamics of the ket and 
the bra and moving only forward 2 ^. Artificial time variables (Keldysh loops) commonly used 
in many-body theor y 24 '? 1 "? 2 , which are necessary for the wavefunction picture, are avoided 
by the density matrix in Liouville space which uses the real laboratory timescale throughout 
the calculation. 
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APPENDIX A: SUPEROPERATOR ALGEBRA FOR THE CURRENT 

A Liouville space operator A a is labelled by a Greek subscript where a — L, R, +, — . It 
is defined by its action on an ordinary (Hilbert space) operator X—. 

(A a X)ij = y^^Ag^a^eX^ (Al) 

A a is thus a tetradic operator with four indices. A^X = AX denotes action from the left, 
and ArX = XA denotes action from the right. We then obtain using Eq. (JAlj) 

(A L )ij, K e = A iK 5ji (A2) 

{AR)ij !K e = A ej 5 iK (A3) 

Note that the order of the j£ indices in Eq. flA3|) has been reversed. 

We further define the symmetric and the antisymmetric combinations A + = \{Al + Ar) 
and A_ = A^—Ar . We then have the tetrahedral matrix elements (A_)ij }K t = A iK 5j£—A£jSi K , 
(A + )ij iK £ = ^[A iK 5je + AijSi K ]. Note that [Al, Br\ = 0. This commutativity of left and right 
operators made possible thanks to the larger size of Liouville space, simplifies algebraic 
manipulations resulting in many useful relations^. 

Using the above definitions of the Liouville space algebra, we can represent / in Eq. (I18|) 
in the form 

1 = 7[ E (T^(*)^(*)(r(^^ 

(r(x 4 )^(x 4 ))„). (A4) 

Here, for a = 1 ip a (xi) = ij}^(xi), and, for a = ip a (xi) — ip(xi). In the r.h.s. of Eq. (jA4|) 
there are sixteen terms. 

We can now expand Eq. (|A4|) in "left" and "right" Liouville space operators^. We note 
that the superoperator corresponding to any function / of an operator Qj can be expressed 
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in terms of Qj + and Qj~, i.e. 



27 



and 



\f(Q. l )\-^f(Q l L)-I((lH) = , I (Q j+ + lQ j -)-f(Q j+ -lQ j -). (AO) 



2 [/(Q,)] + = f(QjL) + f(Qjn) = f (Q j+ + \Qj-) + f (Q j+ - \Q 3 -) ■ (A6) 
Using these relations, we have 

(V> Q (zi)V> M/ Vi))- = (^i)^ M/ Vi))l - {r{xx)iP m {x x )) R 

= rM^ix,) - v#Vi)vs(*i), (A7) 

substituting Eq. ()A7|) in Eq. ()A4|) . we obtain 



r],K,\,v=L,Ra^!3=0,l 



m^f p {x,)r v {x A )^{x A ))) (as) 

Eq. ()A8|) which has 196 terms is used to derive Eq. (|19|). 

APPENDIX B: RELATIONS BETWEEN LIOUVILLE AND HILBERT SPACE 
GREEN FUNCTIONS 

In this Appendix we establish a one-to-one correspondence between Liouville and Hilbert 
space Green functions. 

Standard field theory is formulated in terms of four types of Hilbert space Green func- 
tion o 24 i 31 i 32 . Using the notation of Ref.— , they are 

G c {x,x) = -i(Ti/>(x)i/>\x)) 

= -i6(t-t)(ij(x)^(x)} + i6(i-t)(^(x)^(x)}, (Bl) 

G d {x,x) = -i(fijj(x)^{x)) 

= -i6{t-t)(i){x)il)\x)) +ie{t-f){^{x)il){x)) ) (B2) 

G > (x,x) = -t(ip(x)^(x)}, (B3) 
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G < (x, x) = i(ip^(x)ijj(x)) 



(B4) 



G c , G c , G > and G < are known as the time-ordered, antitime-ordered, "greater" and "lesser" 
Keldysh Green functions, respectively. T is a time-ordering operator in Hilbert space, which 
takes any product of Fermi field operators ( and \ and reorders them in ascending times 
from right to left, i.e. 

C(Tl)x(72) T 2 < 71 

T((n) X (T 2 ) = - X (r 2 )C(r 1 ) n < r 2 (B5) 

. l^Mn) - x (n)an)} r 2 = n 

T is an antitime-ordering operator in Hilbert space, which reorders any product of Fermi 
field operators in descending times from right to left, i.e. 



f C(n)x(r 2 ) = 



-x(t2)C(ji) r 2 < n 

C(n)x(n) ti < t 2 

{ |[C( T i)x( T i) — x( T i)C( T i)] r 2 = r 1 



(B6) 



Four types of the Green functions naturally show up in Liouville space: 

G LL (x,x) = -i{T{^ L {x)^\{x)))\ 



(B7) 



G RR (x,x) = -i(T(ip R (x)ip R (x))}; 



G LR (x,x) = -i(T(ip L (x)il) R {x))); 



(B8) 
(B9) 



G RL (x,x) = -i(T(iJ) R (x)iJ)l(x))), 



(BIO) 



where T is a time-ordering operator in Liouville space, which rearranges all superoperators 
so that time decreases from left to right, i.e.— 

CA t i)x^(t 2 ) t 2 < n 

TUn)x,{r 2 ) = <j -x^Un) Tl < r 2 (Bii) 

t - Xn(n)Cu(n)] t 2 = t x 

where Cuij) is a Liouville space Fermi superoperator, and v and fi can be either +, — , or L 
or R. 
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We next establish the connection between the Hilbert space (Eqs. f)Bl|) - (JB4J)) and Liou- 
ville space (Eqs. (jB7j) - (|B10jl ) Green functions. For Grl and Glr we have 



G RL (x,x) = -iTr 



fi/> R (x)ipl(x)p 



-iTr 



—iTr 



ip(x)i) ] (x)p 



G > (x, x), 



(B12) 



G L r(x,x) = -iTr Tip L (x)ip R (x)p 
= G < (x, x), 



-iTr 



tp(x)p ip i '(x) 



-iTr 



^(x)ip(x)p 



(B13) 



where po is unperturbed many-electron density matrix of the isolated molecule uncoupled 
to the electrodes. 

We next turn to G RR and Gll- For x > x, we have 



Grr(x, X) 



-iTr 
-iTr 



fijj R (x)^ R (x)po 
i)\x)i)(x)po 



-iTr 



p ^{x)i/;(x) 



(B14) 



G LL (x,x) = -iTr Tip L (x)ipl(x)p 



—iTr 



ip(x)^(x)p 



For x < x we analogously obtain 



G RR (x,x) = -iTr ifj(x)^(x)p 



(B15) 



(B16) 



G LL (x,x) = -iTr ip\x)il)(x)po 



(B17) 



Comparing Eqs. ()B14|) . ()B16|) with Eq. ()B2|) . we see that G RR (x,x) = G c (x,x). Comparing 
Eqs. (IB15J) . ijBTTj) with Eq. (pljl . we get Gll(x, x) = G c (x, x). Eqs. (IBTl and (IBTal) provide 
the other relation G R l(x,x) = G > (x,x) and Gl R (x,x) = G < (x,x). 



APPENDIX C: DFT KELDYSH GREEN FUNCTIONS IN HILBERT SPACE 

The four Hilbert space Green functions (Eqs. (|B1|) - ()B4J0 are not independent. To ob- 
tain the relations between them we introduce the retarded G r and the advanced G a Green 
functions^ 

G r {x,x) = 6{t-i){G > (x,x) -G < {x,x)], (CI) 
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G a (x,x) = 8{t-t)[G < {x,x)-G > {x,x)\. (C2) 

The retarded DFT Green function in real space is given by Eq. (3.2) of Ref.— (see 
Refs. 24 ^ 34 ' 36 ) 

CT(x, x) = f *!><M> £ (C3) 

V ' J 27T a W- (CJ Q -CJ F ) +«f V ^ 

where 5 — > 0, aj a is the energy of a Kohn-Sham orbital molecular orbital; and Up is the 
Fermi energy of the molecule. 

We define retarded DFT Green function G r a {uS) of the Kohn-Sham orbital a in the fre- 
quency domain 

G r (x,x) = J^e-^^ a (v)^ a (f)Gl(u;), (C4) 

where 

GTJu) = (uj-cOa + UF + iSy 1 . (C5) 
To connect G r \ G a , G > ,G < , G c and G' c ^- we use Eqs. JOB, (102) together with 

G c (x,x) = -G"*{x,x), (C6) 

G a (x,x) =G r *(x,x), (C7) 

G c (x, x) + G%, x) = G < {x, x) + G > (x, x), (C8) 

Since G r (x, x) is given by Eq. ()G3|) . we can use the system of five equations 
Eqs. dHH), (jH2|) . (fHfijl - (jHSjl to obtain all Keldysh DFT Green functions for the molecule 
alone G a , CP, CP, G c and CP Solving these equations, and switching to the Kohn-Sham 
orbital representation in the frequency domain (see Eqs. ()C3j) - (|C5J) ) we obtain all Keldysh 
DFT Green functions for the molecule 

G c a {ui) = (uj + uop — u a + id) 1 9(luf — u a ) + (u + up — ou a — idy 1 9(uj a — up), (C9) 

G„(u) = — (u; + cuf — — 9(uf — uj a ) — {uj + ujf — uJ a + i5)~ l 

9(uj a -uj F ), (CIO) 
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G a a (oo) = {u - oo a + oo F - td)' 1 , (Cll) 

G>(u) = -2m(l-0(u a -u F ))6(w-u a + u F ), (C12) 

and 

G^(u) = 2m9(uj a — ujf)8{u — u a + top). (C13) 

The Keldysh Green functions for the electrode A can be obtained analogously by replacing 
G r (x,x) with D r A (x,x). The retarded Green function of an electron D r A (uj) in a metal 
electrode in frequency domain is 

D A (u})~(u-(i A + i5)- 1 , (C14) 

where ^a is the Fermi energy level of electrons in electrode A (B). The difference between the 
chemical potentials of the electrodes A and B is given by the external voltage U = fiA~ Hb- 
In Eq. (|C14J) we assumed that the electrons in the metal are at the Fermi level. 

D A (u) = {lo - \la + iSy 1 9(uj - //a) + - AU — *<5)~ 1 

e(fiA-oo), (C15) 

D A (u) = — (oj — Ha — iSy 1 9{u - fiA) — (w — ha + 

9{va-u), (C16) 

D%(u) = (u-n A -i6)-\ (C17) 

D>(uj) = -2m(l-6(u-n A )), (C18) 

D%(u) = 2ni6(u-n A ). (C19) 

Similar expressions apply to electrode B. The final expression for the molecular current 
only depends on the sum of the Green functions of electrodes A and B (D c = D C A + D C B , 
D' C = D\ + D%, D> = D> + D> and D< = D< + D<). 
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APPENDIX D: THE VARIOUS CONTRIBUTIONS TO THE MOLECULAR 
CURRENT 

Only terms which contain the contractions (T(^) L (x n )il;' L (x m ))), {T{ipR{x n )ip R {x rn ))) , 
(T(ijjR(x n )ijj] L (x m ))) , an d (T(tfj L (x n )ifj R (x m ))) (and similar contractions for the electrodes) 
can contribute to / in Eq. (fTT?|) . because all other contractions vanish. Therefore, only 
twenty one terms survive in the r.h.s. of Eq. (|TT?|) . These terms contain the four types 
of Liouville space Green functions of the molecule Gll, Grr, Gl R , and Grl introduced 
in Appendix El and of the electrodes D LL = Dll(A) + Dll(b), D rr = D RR(A) + D RR(B) , 
Dlr = D LR ( A ) + -Dlr(b), and D RL = Drl(a) + D RL ^ B y 

As an illustration we write the expression for one of the terms contributing to S, 

J15 = -ill (x 1 ,x' 1 )U(x 2 , x' 2 )U(x 3 , x' 3 )U(Xi, x'^Grl^X!, x 2 )G lr (x 3 , x a )G ll (x, x) 

D lr (4^'i)Drl(x^x' 3 ). (Dl) 

Using the relations of Appendix El we can recast this in the form of Eq. (jD16|) . 

Below we present all twenty one terms which contribute to the current (Eq. (|22)l). using 
the Kohn-Sham orbitals representation in the frequency domain. The Greek indices run 
over the occupied Kohn-Sham orbitals a, 0, 7 = 1 ... N. 

h = ^E f/ ^^>(^)G^(a; 1 )G c J^)^(-^) J D<(-^); (D2) 

a 

h = ^E^^G'>(c)G^( Wl )G^(^) J D<(-c 2 ) J D c *(-c 1 ); (D3) 

a 

h = -iY,Up a G>(u)G<(u l )G>(u 2 )D^(-u J2 )D < {-uj 1 )- (D4) 

a 

h = zE^^^(^ 2 )G^(a;)G>(cu 1 ) J D c *(-a; 2 ) J D < (-^); (D5) 

a 

h = -^E f/ /3f / '^( Wl )G>MG^(^ 2 ) J D < (-^) J D c *(-^ 1 ); (D6) 

a 
a 

h = -«E^^G<( W )G:( Wl )G:( W2 )r(- W2 )5>(- Wl ); (D8) 
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a 

-^E^f / ^<MG>(^ 1 )G<(^) J D*(-c 2 ) J D>(-c 1 ); 

a 

iE^^^(, 2 )^(,)G<( Wl )r(-a; 2 )i;>(-a; 1 ); 

a 

a 

-iY.UplG^u^G^Gliu^i-u^i-u,). 

a 

^[/^ T ^G;( W )^( Wl )G;( W2 )r(- W2 )D w (- Wl ); 

cry 

-^^C/ 7 ^G'^(u;)G'^(a; 1 )^(a; 2 )L' c *(-a; 1 )L' £ *(-a; 2 ); 
-j^f/^ 7 {/%MG<WG>WD<(- W2 )D>(- Wl ); 
-tJ2U^UlGUuJ2)G<(u 1 )G>(u)D"^-u 2 )D>(-u 1 ) ] 

cry 

^U p U,UlG^u; 1 )Gi(u;)G>MD<(-u; 2 )D>(-u; 1 y, 

cry 
cry 

-^E^f / 7^<(c)G>(cc; 2 )G;(c 1 ) J D<(-c 2 ) J D c *(-c 1 ); 

cry 

-^Ef / ^7^^(c 1 )G>(cu 2 )G<M J D<(-c 2 ) J D £ *(-c 1 ); 

«7 

-i^(7,[/ 7 ^G<HG:MG>( Wl )^(- W2 )Z)<(- Wl ). 

Q?7 
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APPENDIX E: RELATION BETWEEN BOSON GREEN FUNCTIONS IN LI- 
OUVILLE AND HILBERT SPACE 



The relations between Liouville space Green functions and Keldysh Hilbert space Green 
functions for fermions also apply to boson fields. We consider a system of identical atoms in 
an external potential. The definitions for the Liouville space Green functions Gll, Glr, 
Grl and Grr (Eqs. (jB7|) - (|B10|) ) are the same for bosons as for fermions with the only 
difference that a time-ordering operator in Liouville space Tg is defined as^ 



Cu{ri)x^2) r 2 < 7i 

xMUn) n < r 2 (El) 

t j[(v(n)x»(n) + Xn(n)(v(n)] r 2 = n 

where („(t) is a Liouville space Bose superoperator, and v and /i can be either +, — , or L 
or R. 

The four types of Hilbert space Green functions for bosons are defined as^ 



G c (x 1 ,x 2 ) = -i(T B 7p(x 1 )^(x 2 )) 

= -ie(t x - t 2 )(ip(x 1 )^(x)) - %6{t 2 - t 1 )(^ t (x 2 )V(xi)), 



(E2) 



G d (x,x) = -i{f B iP(x)^(x)) 

= -i9(t 6 - t 5 )(tfj(x)^(x)} - i6(t 5 - t 6 )(^(x)tfj(x)} 7 



(E3) 



G > (x 1 ,x 2 ) 



-i{^(xi)^(x 2 )), 



(E4) 



G < (x 1 ,x 2 ) 



-i(ip ] \x 2 )ip{x 1 )) . 



(E5) 



In Eq. (|E2|) Tg is a time-ordering operator in Hilbert space for bosons, which takes any 
product of Bose field operators (, x and reorders them in ascending times from right to left, 
i.e. 

C(n)x(r 2 ) r 2 < n 

T B ari)x(r 2 ) = X (r 2 )C(r 1 ) n < r 2 (E6) 

J[C(nMri)+x(ri)C(ri)] r 2 = n 
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Tb is an antitime-ordering operator in Hilbert space for bosons, which reorders any prod- 
uct of Bose field operators in descending times from right to left, i.e. 

xfaXfa) r 2 < n 

T B C(ri)x(r 2 ) = C(n)x(r 2 ) n < r 2 (E7) 

J[C(ri)x(ri)+x(ri)C(ri)] r 2 = n 
In analogy with Appendix FBI we can show that for bosons Gll{x\, x 2 ) = G c (xi,x 2 ), 
G RR (x 1 ,x 2 ) = G c (x 1 ,x 2 ), G RL (x 1 ,x 2 ) = G > (x 1 ,x 2 ), G LR {x x ,x 2 ) = G < (x 1 ,x 2 ). 
The retarded boson Green function in real space is given by22 

crfr,*,) = /^e'^-. E ^("WM (Eg, 

where 5 — > 0; w„ is the energy of a eigenvalue; fi is the chemical potential; y a ( r i) is the 
wavefunction corresponding to a eigenvalue of the external potential. 

We define retarded boson Green function G r a {uS) in the frequency domain 

G r (x u x 2 ) = /^e-^^)E^( r i)<( r 2 )^M, (E9) 

where 

G r a {u) = (u-Ua + v + iSy 1 . (E10) 
The Keldysh Green functions for bosons are 32 

G c a {uj) = (u + // — u a + idy 1 — 27rin a 5(uj — uj a + /i), (Ell) 

G c a {u) = — (w + fi — uj a — id) 1 — 2mn a S (a? — u a + /i), (E12) 

G a Ju) = (uj-u a + fi)- i5Y l , (E13) 

G>(lo) = -2Tii(l + n a )5(uj-uj a + ij), (E14) 

G^(w) = — 2TTin a 5(uj — uj a + /x), (E15) 

where n a = (exp[(u; Q , — /i)/(/ceT)] — 1) _1 is the occupation number of the a harmonic 
oscillator (T is temperature, and k R is Boltzmann constant). 
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